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We present an algorithm for the derivation of Dyson-Schwinger equations of general theories that is 
suitable for an implementation within a symbolic programming language. Moreover, we introduce 
the Mathematica package DoDS$\ which provides such an implementation. It derives the Dyson- 
Schwinger equations graphically once the interactions of the theory are specified. A few examples 
for the application of both the algorithm and the DoDSE package are provided. 
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I. INTRODUCTION 



Correlation functions are the basic quantities in local quantum field theories and encode all physical information 
about the theory. They fulfill quantum equations of motion, conventionally called Dyson-Schwinger equations (DSEs) 
^ ' PL 0] which are related among each other and form a set of infinitely many coupled equations. Derived from the 
translational invariance of the path integral they are genuinely non-perturbative and describe the physics of the system 
on all scales. This makes them a very useful tool for investigating aspects on which some alternative approaches fail. 
Perhaps the most prominent example is perturbation theory, which is not valid in the strong coupling regime. Since 
DSEs are likewise applicable in the weak coupling region they successfully extend the results of perturbation theory 
into the strong coupling domain. An alternative non-perturbative tool, which can be used complementary to DSEs, 
are Monte-Carlo lattice simulations. Due to the discretization of space-time they have their limits for very low and 
very high momenta, the former being restricted by the size of the lattice and the latter by finite lattice spacings. 
DSEs, in contrast, are formulated in continuous spacetime and allow to study also the analytic structure and the 
infrared regime which is particularly important in an asymptotically free but confining gauge theory like quantum 
chromodynamics (QCD). 
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However, DSEs also have their challenges. They represent strongly non-linear integral equations that are numerically 
involved. Moreover, as they form an infinite tower of equations, they have to be truncated. Recently it turned out 
in the context of Landau gauge QCD that the leading order truncation based only on the propagator DSEs can miss 
important qualitative features that are encoded in the equations for the vertices. In particular, the quark-gluon vertex 
provides a novel mechanism for confinement and chiral symmetry breaking [jj, as well as anomalous mass generation 
[4J in QCD. Yet, the DSEs for the vertices become increasingly complicated and correspondingly hard to obtain 
algebraically. Another complication is given by the necessity of gauge fixing and the additional degrees of freedom 
and interactions arising from the corresponding constraints. In particular, in non-covariant gauges, like Coulomb 
gauge or non-linear gauges like the maximally Abelian gauge, this increases the effort to derive the DSEs already at 
the propagator level considerably, see for example 0, S| , and calls for an algorithmic way to derive these fundamental 
equations. This is especially useful when working with actions that contain many different fields and interactions, 
as arise e.g. when symmetries are not manifestly realized or in the case of necessary additional terms in the action. 
Examples for the latter are the four-ghost interaction required to ensure renormalizability in maximally Abelian gauge 
2, B or generalized constructions of Lagrangians allowing additional terms as in ghost-antighost symmetric gauges 
g.flQl|. A convenient way to derive DSEs also simplifies the comparison of different gauges required to obtain a more 
gauge independent picture of the basic underlying mechanisms. Finally, it is particularly useful in the context of an 
IR analysis where the IR scaling, which is important for long-range properties like confinement, can be abstracted 
from mere power counting. 

The aim of this paper is to present such an algorithmic derivation of Dyson-Schwinger equations. A similar aim 
has been extensively followed in perturbation theory where it resulted in the basically automatic computation of 
numerous physical processes to a given order, cf. e.g. [TJ [12, US, El, EH]. Here we partially extend this idea to 
the non-perturbative regime were such an automatic solution of the created equations is surely beyond our scope. 
Instead we present an algorithm for the derivation of the equations that is suited for implementation into a symbolic 
programming language. This algorithm is presented below in Sec. lUand is implemented in the Mathematica package 
DoDSE, which stands for Derivation of Dyson-Schwinger Equations. An example of how to use the algorithm is 
provided in Sec. IIIII Details on the DoDSE package are presented in Sec. IIVI Whereas a direct algebraic derivation 
can be quite a tedious task, with the symbolic and graphical notations employed, one can obtain DSEs for general 
actions with a relatively high number of interactions. Moreover, the presented algorithm operates directly on the 
level of the effective action and circumvents the tedious step to decompose connected into proper vertices necessary 
in a derivation on the level of the generating functional of full Green functions. We implemented this algorithm 
up to the diagrammatic level into DoDSE. From the interactions of the theory, given in symbolic form, the code 
derives DSEs to the desired order. The outcome are symbolic representations of the Feynman diagrams encoding 
their topological structure and their symmetry factors. The last step, which has to be done manually, to get the 
full algebraic form of the DSEs is the replacement of the symbolic form by the explicit integral expression involving 
proper and bare correlation functions. For some applications of DSEs it is not even necessary to process the symbolic 
equations even further since they can be used directly, as is e.g. the case for scaling analyses. Finally, we note that 
the presented algorithm is in principle also applicable for the generation of a perturbative expansion by re-inserting 
the Dyson-Schwinger equations for dressed vertices and truncating at a given loop order. 



II. DERIVING DYSON-SCHWINGER EQUATIONS IN SYMBOLIC NOTATION 

The algorithm presented here involves several abstractions. The first is the use of a symbolic notation that enables 
us to keep expressions relatively short. Writing out the equations in full detail, i.e. Lorentz and internal group indices 
as well as all coordinates in position space or momenta, leads to expressions that can easily hide the underlying 
basic structure. Secondly, we employ a super-field formalism that includes all irreducible fields of the theory into a 
single reducible multiplet. Finally, we exploit the fact that a local action can be expanded in the fields. Thereby the 
expansion coefficients, that can be operators in coordinate as well as internal space, can be left unspecified during 
the derivation and only have to be inserted in the end to obtain algebraic expressions for the graphical equations. 
Thereby, the main steps can be done using diagrammatic replacement rules. All the techniques used in the following 
are not new, but combined they provide a powerful way of deriving DSEs conveniently and fast. As an explicit 
example we demonstrate the procedure for the two inequivalent DSEs of the ghost-gluon vertex in Landau gauge 
Yang-Mills theory in Sec. IIIII Before we explain the actual algorithm in Subsection III Bl we make a short summary 
of used identities and employed conventions. 
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A. Basic Identities and Conventions 



The basic object we will start from is the generating functional of one-particle irreducible (1PI) Green functions. 
In the following we do not specify any particular fields but will rather use a more compact notation with a super- 
field denoted by fc. It represents a reducible multiplet containing all irreducible representations present in the 
action. We also introduce corresponding sources J{. The multi- index i includes the label of the particular irreducible 
representation, its internal indices and also the space-time dependence of the field. For example in Landau gauge 
QCD 4>i is the set of the gluon, ghost and quark fields: {A^(x), c a (x), (?(x), q l (x), (?(x)}. Where appropriate the 
Einstein convention for summation is understood and amended by integration. 

The effective action is given by 

r[$] = -w[j] + $iJi , (i) 

where W[J] is the generating functional of connected Green functions related to the functional path integral by 

Z\J\ = f D[4>]e- s+tt >> J > = e w ^ (2) 



with S being the gauge-fixed action. The effective action r[<f>] depends on the averaged fields <& in the presence of 
external currents J, 



$, = (Mj = -fr = ZiJ]- 1 I DWke-s+M* . (3) 



SW 
6Ji 

For Grassmann variables terms like (f>jJj include the fields and sources in the usual order 1 , i.e. for example in Landau 
gauge QCD we have fajj = J dx(A^(x)j^(x) + a a (x)c a {x) + c a {x)a a {x) + ff{x)q i {x) + tf(x)rf(x)). 
For non-vanishing external sources (denoted by the index J) the propagator of the multiplet <j) is 

Al^^=(J^-Y 1 (4) 



but for vanishing external sources only some elements remain. For instance in Landau gauge QCD these would be a 
diagonal element for the gluon and off-diagonal elements for ghosts and quarks, but also mixed propagators can arise 
if allowed by the action. For the derivation of DSEs it is important to keep the general expressions and only at the 
end the external sources can be set to zero. Otherwise one would miss contributions. 
Higher 1PI vertices are defined as derivatives of the effective action: 

r J = — (5) 

In the following we will need the derivatives of fields, propagators and vertices with respect to fields. These are given 
by 

^7$i = Sa , (6a) 



s f s 2 r \ _ / 5 2 r y 1 ( a 3 r \f s 2 r 

5 „ r 5T 



-r J ■ = — = r J - ■ (6c) 

and are represented graphically in Fig. [TJ 

The physical correlation functions are obtained from the above expressions when evaluated at the vacuum expec- 



1 Another possibility would be the use of a metric as for example used in ref. |T3 |. 
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FIG. 1: Diagrammatic rules for differentiating an external field, a propagator or a vertex. The circle with the cross denotes an 
external field, small blobs denote dressed propagators, and big blobs 1PI vertices. The double line represents the super-field </>. 



tation values of the fields, corresponding to vanishing external currents, 



r 



h " ' ( Pin)lPI = 



A J | 



<E>=<I> n 

ST 



= r 



J 



where the vertices involve only proper, i.e. 1PI, diagrams. By construction the arising generating equations always 
involve one bare vertex function 



K ••■</Oi 



ss 



PI 



$ = <f>0 



sw 
~sj1 



(7) 



,7=0 



We comment on the inclusion of fermions in Subsec. Ill CI 



B. Algorithm 



The functional DSEs [171 . lift |l9|| are obtained by the invariance of the path integral with respect to variations of 
the integration variable <j> 2 : 



-Z[J] = J D[<t>] (-^ + J,) e -*+^ = 
+ Ji Z[J] = 0. 



SS 



(8) 



4,'.=S/sji 



Substituting Z[J] by e w ^ and using 



-W[J] 



SJi 



,W[J] 



8W[J] _S_ 
SJi SJi 



after multiplication of eq. ([8|) from the left with e w ^ we find 

SS 



Ji = 0. 



(9) 



(10) 



This is the functional DSE for connected correlation functions. To get the corresponding version for 1PI functions we 
perform the Legendre transformation of W with respect to all sources. Thereby SW[J]/SJi changes to $i and S/SJi 
becomes 



_S__S$ 1 _S_ _S_SWJ_ S 2 W S 5 
SJi~ SJi 5<f>j ~ SJiSJj S<f>j ~ SJiSJjWj ~ lj W 3 



(11) 



2 For a different approach using equal time commutation relations and Heisenberg's equation of motion see ref. |2C| ]. 
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FIG. 2: The functional DSE for 1PI functions. Crosses in circles denote external fields. All internal propagators are 1PI and 
the big blob denotes a 1PI vertex function. 



This yields 

SS 



fo=*i+A^ (5/5*3- 



We stress again that here the summation over the index j includes summation over different fields. To simplify the 
process of differentiating the action we expand it in the fields 3 . To stay as general as possible we again use the 
super-field 4> 

S[4>] =yS rs (j) r 4> s - — ] S r st4>r<t ) s<t)t - -^Srstu^r^sMu- (13) 

but for explicit calculations one can alternatively use an expansion in the fields of the action at this point. The Si... 
are the expansion coefficients and correspond to bare quantities (or the inverse ones in case of the propagators) of 
which for an explicit action usually several vanish. Their signs have been chosen in accordance with eq. (7]). For the 
propagators as well as for all derivative interactions these are operators in coordinate space and may have a non-trivial 
structure in the internal space as well. We now differentiate once with respect to a field <f>i, 

SS 1 1 

— = S is (f> s - —Ststfisfit ~ 7^Sirst4>s4>t4>u, (14) 

o<pi 2! SI 

and replace the field operators according to eq. ifl2|) by 

"< • A :'^r;- (is) 

Performing the derivatives and employing eq. © we get the general generating DSEior 1PI Green functions, depicted 
in Fig. 1 



— =s is $ s - - 

1 



S ls $ s -^ st ($ s $ t + A s J t ) 



3 , + + A^A^,.A;[ r l ;'„.,). (16) 

from which the DSEs of arbitrary Greens functions can be obtained. This equation can be decomposed into generating 
equations for the individual irreducible fields which by construction always involve one bare vertex function. If we had 
calculated these equations by direct differentiation of the action as given in eq. (fT2|) , it would have been a tedious task 
to get the bare vertices and propagators in the final expression. In particular, in this case derivatives on (^-functions 
occur that have to be properly resolved by partial integration. This is the advantage the formal expansion in fields 
entails. In this formalism everything is hidden in a single index of the fields and the details of the bare vertices only 
have to be specified in the end in order to obtain algebraic expressions for the DSEs. We have checked by an explicit 
computation in the case of non-Abelian gauge theory that both ways to compute eq. (fIBf indeed reproduce the same 



3 We restrict ourselves here to quartic interactions since they are renormalizable in four dimensions, but there is no restriction on the 
order of the interactions, like e.g. five-point vertices, in the DoDSE package. 
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FIG. 3: The DSE for a generic two-point function. 



result. 

The corresponding DSEs for arbitrary correlation functions are obtained by further functional derivatives of the 
generating equation eq. lfl6|) which are computed via eqs. $6§. From now on we can proceed by use of the diagrammatic 
rules given in Fig. [IJ where all internal lines denote dressed propagators. Here, a major advantage of the diagrammatic 
rules is that one does not have to take care of indices. We would like to stress again that the appearing super- 
propagators have off-diagonal components corresponding to mixing of the irreducible fields. These mixed propagators 
are important in the derivation process, although they seem like an artificial complication. In "simple" cases as for 
example the three-gluon vertex DSE in Landau gauge the result would not change if we kept only "real" propagators. 
However, for the ghost-gluon vertex (see Subsec. IIIII for details), for higher vertex functions and in certain gauges 
already for the propagators, some contributions would be missing. 

The DSE for a generic two-point function is derived by performing another differentiation of the generating DSE 
eq. (fT6|l using the diagrammatic replacement rules of Fig. Q] in the corresponding diagrammatic representation Fig. 
[2] in all possible ways. The result is shown in Fig. 03 Proceeding to higher vertex functions the number of diagrams 
increases rapidly: For three-point vertices there are 15 generic diagrams and for four-point functions 60. For real 
applications it is therefore recommendable to exploit possible simplifications. First, the final number and form of 
graphs depend on the first differentiation in eq. lfT6|) as the corresponding field determines which bare vertices appear 
in the diagrams. For example the DSE of the ghost-gluon vertex in Landau gauge QCD has only four terms, when 
the first derivative is performed with respect to a ghost field. In this case one can drop all diagrams with bare gluonic 
vertices. On the other hand, if one starts with the gluon field, all vertices have to be kept and one ends up with twelve 
graphs. Secondly, one can skip some diagrams taking into account where one is going. Simple examples are that for 
a three-point function we do not have to drag along the bare four-point vertex or diagrams with an external field can 
be dropped if no further derivatives with respect to this particular field follow. 



C. Inclusion of Grassmann Fields 



Anticommuting fields need slightly more care when performing derivatives to get the correct signs. In this subsection 
we denote Grassmann fields by ip and ip with sources fj and rj, respectively. Let us first specify our convention for the 
derivatives with respect to Grassmann quantities. We choose left- and right-derivatives: 

A •= 1 ± •= 1 ( 17) 

5ip ' 5tp' Sip 6ip K ' 

This entails that the proper definition of a Grassmann field propagator is 

$ 2 W ( S 2 T V 1 

A« = -==-=- = 7TT- ■ (18) 



13 SrjiSrjj \Sipi6ip 

In general derivatives always have to be ordered such that derivatives with respect to Grassmann fields are right of 
those with respect to anti-Grassmann fields. A quartic Grassmann interaction then has the form 



I ijkl — 



dipidipjdipkdipi 



(19) 



For easier readability and also in correspondence with DoDSE the indices of the T do not reflect the order of how the 
derivatives are performed, but rather have the order in which the derivatives appear, i.e. that the differentiation with 
respect to ipj has to be performed before that with respect to ipi, but ipk comes before ipi. 
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= T J - ■ 
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-r J 
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= T J 
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When using the algorithm described above derivatives always act from the corresponding direction. Eqs. J6]) have 
to be amended by 

ik, (20a) 

it, (20b) 

^X^Ji-Jn = (20c) 

where here the superscript J denotes the dependence on all sources and the additional superscripts of the r J denote 
the fields corresponding to the indices. This means in turn that at the end, when the external sources are set to zero 
and the reducible multiplet is decomposed into the irreducible, physical fields, some Grassmann derivatives can be 
unordered. Ordering them gives the signs expected normally for Feynman diagrams with fermion loops. However, 
in some cases this algorithm is oversimplified when several super-fields are involved. Using diagrammatic rules this 
problem can be circumvented and the minus signs for closed fermion loops have to be added manually at the end. 
Since DoDSE needs specific rules how to perform the derivatives, it may happen that the wrong sign appears for 
two-loop graphs. 

A simple example for the change of sign is the quark respectively ghost loop in the gluon DSE of Landau gauge 
where tp = {q, c} . Replacing ip r by f/v + in the corresponding part of the first derivative of the action with 

respect to Ai, 



S 



ITS 



VVV> S -> -Stf" (VvV's + A r J ^') , (21) 



and differentiating once more with respect to Aj yields 



-S£t>A++A+? - - f T - - + gluonic terms (22) 



Ordering the derivatives changes the sign of the expression and leads to the expected relative minus sign of closed 
fermion loops. 

Finally we should comment on the expansion of the action when Grassmann fields are involved. First, we recommend 
to order Grassmann interactions such that all anti-Grassmann fields are left of the Grassmann fields. Second, the 
expansion coefficients are antisymmetric in the indices belonging to Grassmann fields. This entails that we can 
differentiate with respect to Grassmann fields as usual, e.g. 

±Sff t m r ^ u = S^HM U ~ StUHr^u = + Stt^rMu = 2Sf^HM u . (23) 

dtpi 



III. AN EXAMPLE: THE GHOST-GLUON VERTEX IN LANDAU GAUGE YANG-MILLS THEORY 

For the derivation of the ghost-gluon vertex in Landau gauge (2ll | we start with Fig. [3l As mentioned above, 
there are two DSEs for this vertex. These arise since the DSEs are derived from the invariance of the initial path 
integral under changes of the fields. The invariance for the individual fields in the theory yields different generating 
equations from eq. fT6j) that can by appropriate functional differentiation generate topologically distinct DSEs 
for the same mixed correlation function. Without approximations these equations should provide identical results 
but they may be affected differently by truncations. Moreover, the dynamics can be represented differently in the 
various equations as demonstrated for the quark- gluon vertex in 0]. We start here with the simpler equation for the 
ghost-gluon vertex, i. e. we take the derivative in eq. (fT6| with respect to the ghost field. Fig. [3]reduces in this case to 
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The symmetry factor 1/2 vanished because there were two possibilities to plug the bare ghost-gluon vertex into the 
third diagram. This is a nice example how the computation simplifies when the possible interactions are reduced by 
symmetries. In the present case the general diagrammatic rules specialize to 



5 A, 



5 A 




which leads to the three-point expression 



| 

p 



| 

= I 




When we set the external sources to zero, the mixed propagators become irreducible gluon and ghost propagators. The 
pure super-field propagator in the second and third terms on the right-hand side yields a sum of different terms, when 
decomposed, but ghost number conservation allows here only vertices with the same number of ghost and anti-ghost 
legs. Therefore, for each diagram only one propagator can be realized. The final result is then the ghost-gluon vertex 
DSE: 



I 



| 
1 

= I 



k 



k 



f X , 

i \ 4 



. cyV2> 



We obtain a second distinct version of the ghost-gluon DSE if we start differentiating with respect to the gluon 
field. In this case the super-field is important as will become evident below. The gluon is involved in all possible 
interactions of Landau gauge Yang-Mills theory, so Fig. [3] is not topologically simplified in this case. For brevity we 
skip diagrams that do not contribute to the ghost-gluon vertex (the tadpole and all incompatible tree graphs as well 
as the graph with the bare four-gluon vertex connected to an external field) . The factor in front of the loop containing 
ghost-fields is changed from 1/2 to 1, because there are two possibilities to insert the bare ghost-gluon vertex as we 
have to consider the direction of fermion lines explicitly. The diagrams left are 



-1 
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Differentiation with respect to the anti-ghost field yields 

k k b b b b b 

ft = I? + i ^ +1/2 ? ^ + / v + /\ + i 1 



✓ v y \ y \ y \ 




+ 1/2 



I 



1/3! 




Again propagators partly involving the super-field are determined by the second field and pure super-field propagators 
by the symmetries of the vertices. Finally the second version of the ghost-gluon vertex DSE is obtained by setting 
the external sources to zero: 



1 



b b b 

+ +1/2 | a + / \ 



I 



+ 



| | | | | | 

+1/2 +1/2 +l/2^# f +1/2 ( <H +1/2 J> } Jfe +1/3 !3' | 1 

% %^ i^-* w 

\ / 

The derivation of the DSEs for the gluonic correlation functions can be done in a similar way. However, since the 
super-field cannot lead to internal ghost loops one can circumvent the super-field formalism and proceed with pure 
ghost and gluon propagators from the level of two-point functions on. For mixed Green functions this would omit 
some diagrams as was explicitly shown for the ghost-gluon vertex, where all diagrams with internal ghost lines - except 
the pure ghost loop (the fifth diagram on the right-hand side) - and the triangle diagrams would be missing. 



IV. DERIVATION OF DYSON-SCHWINGER EQUATIONS USING MATHEMATICA 

As should have become clear the manual application of the algorithm becomes rather tedious for more complicated 
Green functions. The algorithm explained in the Sec. [i] is perfectly suitable for an implementation into a symbolic 
programming language like Mathematica [12] • The functionality of the package DoDSE is that the user enters the 
interactions of the theory and then all desired DSEs are automatically generated from it. In special cases it may in 
addition be necessary to provide further information about symmetries of vertices. Together with the package there 
is also a notebook (examples. nb) available that contains the model used below as well as other examples like Landau 
gauge QCD. 
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A. Technical Notes 

To load the package file DoDSE.m, one either copies it to the subdirectory DoDSE of $UserAdd0nsDirectory 4 
and evaluates « DoDSE 'DoDSE' (which would be the standard way of installation) or uses the Get command to 
load it from any other place: Get [pathToTheFile] . DoDSE was developed under Mathematica 6 and will not work 
with Mathematica 5.2 or lower since some functions new to Mathematica 6 are used. There exists no dedicated 
documentation within Mathematica's Documentation Center, but help on single commands is available using the 
command ?, e.g. ?DSEPlot. 

B. Using DoDSE 

For the derivation of DSEs via the DoDSE package it is usually only necessary to specify the interactions of the 
given theory. This is done by a list whose elements are given by the individual propagators and vertices in the bare 
Lagrangian of the theory. These are in turn represented by lists containing the external fields of the corresponding 
correlation functions. With this information the function doDSE can now derive any DSE. The output can be used 
directly, put into a more convenient form using short notations for propagators and vertices or plotted with DSEPlot. 
The complete calculation takes seconds, maybe minutes for vertex functions with more than four legs. As the package 
is written for Mathematica we have a wide range of tools available for processing the results further, for instance the 
uniform IR scaling exponent of a graph (cf. e.g. [H, 0|) can be derived directly from the output of DoDSE. 

In order to demonstrate the method to derive DSEs let us consider an example model. Its action consists of four 
fields, two of them bosonic (A, B) and two fermionic (c, d): 

C = — Sij AiAj + 2^tj BiBj + Sfj CiCj + Sfj d~idj-\- 

- S:)yA,c,n, - ^S^AiAjBk \s^ BB AA^M ^S^A^A.Ai - S^fcdjdm . (24) 

We wrote the quartic Grassmann interaction such that no additional minus sign occurs when differentiating with 
respect to c and c or d and d. The interactions of this action are entered in the form 

ilist = {{A, A}, {B,B>, {cb,c}, {db,d}, {A,cb,c>, {A,A,B>, {A,A,B,B}, {A, A, A, A}, {cb,db,d,c}}; 

With this representation of the theory, we can start deriving the propagator DSEs using the function doDSE which 
generates the non-trivial right hand side of the DSE for the corresponding correlation function. As arguments it takes 
the list of interactions ilist, and a list of the fields with their respective indices included in the correlator for which 
we want to derive the DSE. In general, the order of the elements in the list of the correlation function for which the 
DSE is derived determines the order in which the individual functional derivatives are taken. When different fields 
are involved this can result in distinct DSEs, as was the case for the ghost-gluon vertex discussed in the last section. 
The corresponding commands for the propagators read: 

AADSE = doDSE [ilist, {A, A}]; 

BBDSE = doDSE [ilist, {B, B>] ; 

ccDSE = doDSE [ilist, {{c, i>, {cb, j}}] ; 

ddDSE = doDSE[ilist, {{d, i>, {db, j}}] ; 

It is not necessary to give the indices of the fields, but it can be done as shown above for the two fermion DSEs. 
Alternatively one can derive several DSEs by passing a list to doDSE, i.e. the following line derives one DSE for all 
primitively divergent vertex functions: 

DSEs = doDSE [ilist, ilist]; 

We can bring the rather long output into a more readable form with shortExpression: 

shortExpression [AADSE] 

or equivalently 



4 Under Unix systems this is normally /.Mathematica/ Applications 
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sE [AADSE] 

oAA_ 1 /pAAAAaAA \ _ 1 /oAABBaBB \ _ cA A BpA ABaAAaBB , qA cb c r A cb c a c cb a c cb 
ij 2^ ij rlsl rl slj 2^ ijrlsl '"Vl sl) J i rl si 1 j tl ul ^rl tl ^sl ul °i rl si 1 j tl ul ^sl tl ^ul rl 

1 /qAAA A r A AAA a A A a A A a A A \ _ ^ j qA A B BpA ABB a A A aB B aB B \ 
g \°i rl r2 si 1 j s2 t2 u2^rl s2^r2 t2^sl u2j ^ V » rl r2 sl J s2 t2 u 2 rl s2^r2 t2^sl u2j 

? AABB r ABA r AAB a A A a B B aB B A A A cA A B BpA B B pA B B a A A a B B aB B aBB 

°\ rl r2 sl 1 j u2 vl 1 s2 t2 ul^rl s2^r2 u2 AA sl ul^t2 vl °\ rl r2 sl 1 j u2 vl 1 s2 t2 ul^rl s2^r2 u2^sl ul^t2 vl 

1 ( oA A A ApA A ApA A A a A A a A A a A A a A A \ _ _ ( qA A B BpA A A r A B B a A A aB B aB B a A A \ 

2 \°i rl r2 sl 1 j s2 tl 1 ul v2 wl^rl s2^r2 v2 AA sl wl^ul tl j ^ V > rl r2 sl J s2 rl ul v2 wl rl s2 r2 v2^sl wl^ul tl j 

1 ( oA A A ApA A BpB A A a A A a A A a A A aB B \ _ }_ I cA A B BpA A BpB B B a A A aB B aB B aB B \ 
r, V°i rl r2 sl 1 j s2 tl 1 ul v2 wl^rl s2^r2 v2 AA sl wl^ul tlj 2 V i rl r2 sl 1 j s2 tl 1 ul v2 wl^rl s2^r2 v2 AA sl wl^ul tlj 

(25) 

shortExpression uses the Mathematica function Style and accepts its corresponding options, e. g. colors 
or FontSize. The symbols used for propagators and vertices are set with the variables $bareVertexSymbol, 
$vertexSymbol and $propagatorSymbol. The standard settings are 5", T and A. The subscripts of these expressions 
are the indices, whereas the corresponding fields can be found in the superscript. Alternatively one can plot the DSEs 
with DSEPlot. To improve the representation it is advantageous to define a few graphics primitives for each field in 
a list of the form {{f ieldl , primitivesl]- , {field2, primitives2} , ...}: 

fieldRules = {{A, Red}, {B, Green, Dashed}, {c, Blue, Dotted}, 
{d, Orange, Dashing [{0 . 02 , 0.01}]}}; 

These primitives can be used in DSEPlot: 

DSEPlot [AADSE, ilist, fieldRules] 




Since DSEPlot uses the Mathematica function GraphPlot it inherits its corresponding options which are unfortunately 
quite limited. For example GraphPlot may show the external points at different places for different graphs. It is 
possible to provide specific coordinates for the external points, but then some graphs may appear with overlapping 
internal lines. For many internal lines it can also happen that GraphPlot draws lines above each other. For this 
reason DSEPlot should be considered more a tool for showing the results of DoDSE than a means of creating nice 
Feynman diagrams. However, the latter is also possible to a limited extend using the available options. In the 
generated diagrams all propagators are dressed. 1PI vertices are denoted by a large blob, whereas bare vertices are 
drawn without a blob. Fermionic fields have an arrow denoting their direction. One can invoke DSEPlot also without 
the additional graphics primitives: 

DSEPlot [AADSE, ilist] 

As a result text labels are attached to the individual propagators denoting their type. 

Alternatively one can get a list of all graphs by setting the option output to List. Other options are indexStyle 
and f actorStyle, which determine the styles of the indices and prefactors. 
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C. Implementation 

In the following we will describe the individual operation steps of the function doDSE for the interested reader 
who wants to do a step-by-step calculation for example to check own calculations. Furthermore using the individual 
routines one has several additional possibilities on which we comment below. We also explain the representation of 
integrals in DoDSE. 

The function doDSE performs the following operations: 

1. It converts the list of interactions to the internal representation of the action with generateAction. Alternatively 
one can directly give such an expression as argument to doDSE. 

2. The first derivative is performed directly on the action using deriv, see eq. fl~2|) . 

3. Then the replacement according to eq. fl5|) is done with replacementCalc. 

4. At this point some graphs may appear several times, identif yGraph adds them up to avoid redundant calcu- 
lations and get the symmetry factors right. 

5. The necessary number of further derivatives can be worked out with deriv. 

6. The external sources are set to zero using setSourcesZero. 

7. For Grassmann fields it may be necessary to order the derivatives with orderFermions such that anti-fields are 
left from fields thereby possibly changing the sign. 

8. Due to our definition of the vertex functions the left hand side contains a minus sign so that the final expression 
is multiplied by —1 for all vertex DSEs. 

In the package several objects are introduced for representing Feynman diagrams. The smallest units are fields 
and indices. Typically they are combined in a list like {A,i}, where A is a field and i its index. They can stand for 
external fields or internal ones. These lists can be grouped to form propagators and vertices, denoted by P and V 
(1PI) or S (bare). Finally external fields, propagators and vertices are combined in an object called op, the biggest 
object, which represents individual diagrams. In it indices occurring twice are summed and integrated over. Note 
that a field is external if it is an argument of op, while fields in vertices represent their legs. An example for a graph 
containing a propagator and an external field {A , j } is 

op[S[{A,i},{A,j},{A,k},{A,l}] , P[{A,k},{A,l}] , {A,j}]. 

The index i only appears once and is therefore an external leg. The op-function is the main object used internally for 
the computations and is also returned in the output of doDSE. It can be used like one would expect from something 
representing a graph, for instance one can sum up several of them or it vanishes if one of its arguments is zero. 
Furthermore it splits up if one argument is a sum. For Grassmann fields there is a peculiarity in the notation: To make 
it as easy as possible to read op-expressions, when two indices appear that are summed over also the corresponding 
fields are written identical, although in the notation above there should be the anti-field in the propagator for 
Grassmann fields. This means that the standard Grassmann propagator is defined as 

P[{c,i}, {cb,j}] , 

where c is a Grassmann field and cb the corresponding anti-field. Vertices are defined as expected, i.e. anti-field 
derivatives are left from field derivatives. As inverse bare propagators are denoted by S, the function of bare vertices, 
they have the convention of vertices in contrast to dressed propagators. 

D. Advanced Options 

In some cases the features of DoDSE described so far are not sufficient. For instance the occurrence of mixed 
propagators or symmetries not apparent in the generic expansion of the fields pose more intricate problems. This 
happens for example in the maximal Abelian gauge, where three-point functions involving three diagonal gluons are 
not allowed due to the color algebra but by default they are constructed by DoDSE, which does not know about the 
color structure. These vertices have to be removed by appropriate restrictions. 
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To understand the necessity of vertex test functions and explicit definition of propagators as explained below, it 
may be helpful to know about the use of the super-field in DoDSE. During the derivation of a DSE the super-field 5 
can occur in two places: vertices and propagators. For the latter a list of rules of possible replacements is needed. 
Without further input DoDSE uses only those propagators that appear in the Lagrangian. With these rules the fields 
in the propagators and vertices are replaced. However, as some of the vertices are forbidden by the symmetries of the 
action, further tests, called vertex test functions, are needed to assess if a vertex can exist in the theory. 

In general DoDSE performs two standard tests. The first checks for number conservation of individual Grassmann 
field species. It can be disabled with the option doGrassmannTest->False for example to allow mixing. The second 
one makes sure that bosons with a discrete symmetry <j> ~ 4> in the Lagrangian only appear in vertices that respect 
that symmetry. These standard tests can be amended by user-defined test functions, which take a vertex as argument 
and give True or False. To employ such functions they are given as argument to doDSE (or setSourcesZero). Vertex 
test functions allow also to truncate the system of DSEs, e.g. by forbidding vertices with a certain number of legs. We 
demonstrate this with a truncation of the AA DSE from above by defining a test function that only allows three-point 
vertices: 

vertexTest [a_V] := Length@a==3 
Called with this restriction 

AADSETruncated = doDSE [ilist, {A, A}, vertexTest]; 

the number of diagrams reduces from 13 to 11, because the two sunset diagrams vanished. 

In the case of mixed propagators in the Lagrangian it is necessary to provide the additional option 
specif icFieldDef initions to doDSE, e.g. 

AADSEMixedL = doDSE[{{A, A}, {A, B}, {B, B>, -[cb, c>, {A, cb, c>, {A, A, B>, {A, A, A}}, 
{A, A}, specif icFieldDef initions->{A, B, {c, cb}}] ; 

If this option was not given, doDSE would assume that A and B are fermion and anti-fermion. The list {A, B, {c, cb}} 
tells doDSE that A and B are bosons, c is a fermion and cb the corresponding anti-fermion. 

If the theory allows mixed propagators not present in the Lagrangian, one has to provide a list of all possible 
propagators to doDSE, e.g. 

AADSEmixed = doDSE[ilist, {A, A}, {{A, A}, {A,B}, {B,B}, {c,cb}, {d,db}}] ; 

where we allowed a propagator between the A and B field. The number of terms in the generated DSEs increases 
correspondingly. The presence of mixed propagators has also a drawback: The final result may contain some terms 
several times. To add them up one employs the function identif yGraphs with an additional option: 

AADSEMixedLId = identif yGraphs [AADSEMixedL, compareFunction -> compareGraphs2] ; 

Since compareGraphs2 can take quite long, doDSE normally uses compareGraphs, which is not adequate for the 
situation with mixed propagators. The number of terms reduces as can be checked with 

countTerms/@{AADSEMixedL , AADSEMixedLId} 

In principle it is possible to derive a DSE step by step instead of using doDSE. Thereby one can track every single step 
of the calculation corresponding to the algorithm described above and perform further manipulations, e.g. keeping 
some external fields. We refrain here from an explanation how to do this and refer the interested reader to the 
notebook examples. nb, which is available together with the package. Finally there are a few tools available within 
DoDSE that allow to check the syntax of expressions or to obtain information about the fields. They are listed in the 
appendix. 



B Its standard name is <j>, but that can be changed by redefining the variable $dummyField. As long as one does not do step-by-step 
calculations it should never appear in output. 
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V. SUMMARY 



In this article we presented an algorithm to derive Dyson-Schwinger equations in a convenient way, in which one does 
not have to deal with the usual abundance of indices and integrals. This is achieved by graphical rules for performing 
derivatives that allow a quick and straight forward derivation even of higher vertex functions. We used this algorithm 
in the Mathematica package DoDSE that can give DSEs once the interactions of the theory are specified. It proves 
especially useful for theories with many interactions or higher vertex functions in general as the number of terms 
grows considerably. This should in particular help to analyze gauge theories in different gauges in order to obtain a 
more gauge independent picture of the described physics. 
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Main functions 


Command 


Description 


doDSE[ilist, clist] 


Derives the DSE for the correlation function clist for a theory with 
interactions ilist. 


doDSE[ilist, clist [, props, 
vertexTest , opt s] ] 


vertexTest is a function for determining if a vertex respects the 
symmetries of the Lagrangian. props is a list of allowed propagators 
given in the form {{fieldla, fieldlb}, {field2a, field2b}, 
...}■. doDSE accepts the options specif icFieldDef initions and 
sourcesZero (prevents the replacement of super-field propagators 
and vertices when set to False). 


shortExpression [expr , opts] 
sE[expr, opts] 


Rewrites a DoDSE expression into a shorter form using 
$bareVertexSymbol, $vertexSymbol and $propagatorSymbol for 
representation. Options of Style can be given. 


DSEPlot [expr, 

ilist [,f Rules, len, opts]] 


Plots graphs, expr is an expression containing op functions, ilist 
the list of interactions and f Rules a list of options for plotting in- 
dividual fields, len determines how many graphs are shown in one 
line. If f Rules is not given, the lines are named according to the 
fields. Possible options are: output->List, to get the result in list 
form, and indexStyle and f actorStyle to change the style of the 
indices and the prefactors (e.g. font size or color). 
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Functions 


for the individual computation steps 


Command 


Description 


generateAction [ilist [,f list]] 


Generates the action in internal representation from the 
interactions of the theory given in ilist. For mixed 
propagators flist specifies explicitly the type of fields in 
the form {bosonl, boson2, {fermionl, antif ermionl} , 
{fermion2, antif ermion2} , 


der iv [expr , dlists] 


Differentiate expr with respect to the fields in dlists. 


replaceFields [expr] 


Replaces the fields in expr by the corresponding expressions after 
the first differentiation is done to change from full to 1PI Green 
functions. 


identif yGraphs [expr [ , 
compareGraphs->cf unc] ] 


Adds up equivalent graphs in expr. cfunc can be compareGraphs 
(standard) or compareGraphs2, the latter being necessary for mixed 
propagators but taking longer. 


setSourcesZero [expr , flist [, props, Sets the external fields in f list to zero, i.e. only physical propaga- 
vertexTest]] tors and vertices are left. vertexTest is a function for determining 

if a vertex respects the symmetries of the Lagrangian. props is a 
list of allowed propagators given in the form {{f ieldla, f ieldlb} , 
{field2a, field2b}, ...}. 


orderFermions [expr] 


Orders derivatives with respect to Grassmann fields such that the 
anti-fields are left of the fields thereby possibly giving a minus sign, 
expr is an op-function or a sum of those. Bare vertices are not 
affected by the ordering. 
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Functions for checks and tools 


Command 


Description 


count Terms [exprl 


C^onnts the number of terms annearinff in the expression 

V J V } I A. 1 1 LjO LjXXV^ XXL1XXXI_/V^X V/ 1 \j \j X XXXO CAi yj yJ CAlL XXXw XXX Ij XX V^-/v k-'-L V-'OOlUll . 


fieldq [f] 


Determines if expression f is defined as a field. 


bosonQ [f ] 


Determines if expression f is defined as a bosonic field. 


f ermionQ [f ] 


Determines if expression f is defined as a fermionic field. 


antiFermionQ [f ] 


Determines if expression f is defined as an anti-field to a fermionic 
field. 


checkFields [expr] 


Checks if all fields in the expression are defined 


checklndices [expr] 


Ohecks if an index" annears more often than twice 

v . xxv^ v->x\.o xx cx>xx xxx vxl-ji. ca yj ij v^clx o xxx kji. \j V/ i \j v^xx ijxxclxx \j vv iv-l*i 


checkSyntax [expr] 


Checks if expr has the correct syntax, i.e. op functions only contain 
nrnna frators vertices and fields 


checkAction [expr] 


Checks if all indices appear exactly twice, the syntax is ok and all 
fields are defined. 


checkAll [expr] 


Performs a series of checks on expr (checklndices, checkSyntax, 
checkFields). 


def ineFields [f list] 


Defines the fields of the action that are given in f list as single 
entries for bosons and grouped by braces for fermions. 


$vertexSymbol 


Symbol representing a vertex in shortExpression. Standard value: 

r. 


$bareVertexSymbol 


Symbol representing a bare vertex in shortExpression. Standard 
value: S. 


$PropagatorSymbol 


Symbol representing a propagator in shortExpression. Standard 
value: A. 



